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ABSTRACT 


An implicit finite-difference method has been developed for computing 
the flow in the near field of a fuel injector. This work was done as part of 
a broader study of the effects of fuel injector geometry on fuel-air mixing 
and combustion. Detailed numerical results have been obtained for cases of 
laminar and turbulent flow without base injection, corresponding to the 
supersonic base flow problem. These numerical results indicated that the 
method is stable and convergent, and that significant savings in computer 
time can be achieved, compared with explicit methods. 
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1.0 INTRODUCTION 

This research program involved experiments and analysis to relate 
characteristics of fuel injector geometry to the fuel-air mixing and com- 
bustion. The original experimental program has been only partially completed 
at this time, because of numerous mechanical difficulties in the test facility 
and the instrumentation. References 1-3 give details of the experimental 
plan and also describe most of the experimental results and conclusions 
obtained to this point. In developing computational methods for the turbulent 
mixing and reacting flow, it became necessary to model the time-average rate 
of fuel consumption by chemical reaction at each point in the flow. References 
3 and 4 describe the main results of this effort. 

The present report is concerned with the last phase of the research 
program, which was the development of a mmierical method for computing the 
flow properties in the near field of a fuel injector. This effort was a 
natural part of a study of effects of injector geometry on fuel-air mixing 
and combustion, because the principal effects of changes in injector geometry 
occur in the near field of the injector. This near field is of great interest 
for fuel injector design, because of the need for mixing and combustion of 
the reactants in minimum combustor length. Important features of the near 
field are: (1) regions of reverse flow (as at the base of an injector or down- 
stream of transverse fuel jets); (2) turning of the fuel and air streams, with 
significant transverse pressure gradients; (3) shock waves and expansion waves 
(for supersonic air and/or fuel streams). It is not possible to use the boundary- 
layer form of the conservation equations to compute even the qualitative flow 
behavior near the injector. This is because the boundary-layer equations tell 
nothing about the three flow features named in the preceeding. 



2 


In developing the computational method and in calculating test cases, 

attention was directed toward a fuel injector in which the fuel jet emerged 

from the base of the injector. The method is not restricted to this type of 

injection, however. For example, transverse fuel, injection into a supersonic 

stream could also be computed. No test cases have been computed for such a 

flow field as yet, however. When fuel injection is through the injector base, 

the flow field resembles that of the well-known base flow problem, but with the 

added complications of fuel injection and subsequent mixing and combustion. 

For this reason, the numerical method was first tested against the supersonic, 

laminar base flow calculations of Allen and Cheng^, and then against the super- 
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sonic, turbulent base flow measurements of Lewis and Chapkis . Both of these 
cases involve homogeneous, nonreacting flows. They serve to evaluate the 
capability of the numerical method to compute a flow with shocks, reverse flow, 
and a wide variety of boundary conditions. A sketch of this flow field is 
shown in Figure 1. 

The specific objective of the present work was to develop a finite- 
difference method for computing the flow in the near field of a fuel injector. 
The general approach was to solve the time-average conservation equations for 
mass, momentum, energy, species, turbulence energy, and turbulence dissipation, 
subject to appropriate boundary conditions. These equations were written in 
time-dependent form, and the steady-state solution was taken to be the solution 

for large time. This is usually the more convenient approach from a numerical 

7 8 

viewpoint . The turbulence transport model was that of Launder and Spalding , 

modified for compressibility effects . The main feature of this model is that 

the turbulence velocity scale and length scale (macroscales) at each point are 

determined from modeled equations for k and e . These then yield the eddy 

2 

viscosity (proportional to pk /e) and the other turbulent transport coeffi- 
cients (each proportional to eddy viscosity) . 




A comment may be made regarding the choice of a turbulent transport 
model that requires two additional partial differential equations, instead 
of a simpler mixing- length model. The attributes and defects of these 
approaches to turbulent transport modeling have been much discussed in recent 
years. It seems worthwhile to point out that there are some tradeoffs involved. 
For the mixing-length model, initial conditions are simplified, whereas 
initial conditions for k and e are usually not available and must be 
guessed. On the other hand, it is difficult to specify an appropriate spatial 
mixing length distribution unless a great deal is already known about the flow. 
Finally, personal computational experience suggests that the e -equation is 
not very well conditioned, in the sense that a poor choice of initial condi- 
tions or too coarse a grid may lead to difficulties in obtaining a numerical 


solution. 
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2.0 NUMERICAL METHOD 

9 10 

The numerical method used was based on that of Briley and McDonald ’ 

As the method is implicit in time, it offers the possibility of significant 
reductions in computation time to obtain a steady-state solution. In the 
Briley-McDonald method, backward time differencing is applied to the conserva- 
tion equations. Then these equations are linearized, with respect to time 
level, by a Taylor series expansion about the nth time level. (Flow properties 
are known at the nth time level from a previous calculation) . Central 
differences are used to replace the spatial derivatives. The result is a 
system of simultaneous linear algebraic equations for the flow properties at 
the. (n+l)th time level. An alternating direction implicit (ADI) technique is 
used to solve this system of equations. 

In the present work attention has been restricted to two-dimensional 
plane flow. This was done because, considering the complexity of the flow, 
it seemed prudent to first develop the method for the simpler two-dimensional 
case. 

Figure 2 illustrates the grid into which the flow field is subdivided 
and the corresponding nomenclature. The finite-difference forms of the 
conservation equations were derived by applying each of these equations to 
a control volume surrounding each grid point. This control volume is 
identical with the cell of dimensions Ax, Ay that encloses the grid point. 
(See Figure 2). This approach has distinct advantages. First, it tends to 
force the conservation laws to be satisfied macroscopically, not only in the 
limit as Ax, Ay, At go to zero (reference 7, p. 28). Second, it greatly 
aids in formulating boundary conditions, and helps assure that physically 

5 

meaningful boundary conditions are applied . 
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The derivation of the finite-difference forms of the conservation 
equations is given in Appendix A. 

The linearization about time step n can be illustrated by considering 
the finite-difference equation for conservation of mass. From Appendix A, 

Eq. ( A5 ), at point i j : 

ap/3t'“= - 6^(pu) - 6y(pv) (2-1) 

Backward time-differencing is then applied to 3 p/St, and the right-hand 
side is treated implicitly by evaluating it at the (n+1) level. Eq. (2-1) 
then becomes: 

(p^'^^ - p’^). ./At = - 6^(pu)^|^ - (2-2) 

Eq. (2-2) is nonlinear in properties at the (n+1) level. So also are the 
remaining conservation equations after being differenced in the same way as 
conservation of mass. The full set of conservation equations in this finite- 
difference form is then a set of simultaneous nonlinear difference equations 
for flow properties at the (n+1) level. The solution of this set of equations 
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would require some time-consuming iterative scheme . 

To avoid an iterative solution, the right hand side of eq. (2-2) is 
linearized about the n level as follows: 

(p’^^^ - p^)/At = - 6^ j^(pu)’^ + (3pu/3t)^ AtJ 
- [(pv)" + (Bpv/3t)’^ At] 


(2-3) 
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The quantity (5pu/St)^ is expressed in finite-difference form as: 

(Bpu/dt)^ = (p"'*'^ - p’^)u’^/At + p’^Cu"'*'^ - u^)/At (2-4) 

Using this method for (Spv/3t)^ and collecting terms, eq. (2c3) reduces to: 

, n+1 n. , n+1 n . n n+1 n n. 

(p - p )/At =-6(p u+pu -pu) 

c , n+1 n . n n+1 n n^ 

- 6y(p V + p V - p V ) (2-5) 

The subscript ij is to be understood in eqs. (2-3) to (2-5). 

✓ 

Equation (2-5) is linear in (n+1) "level variables. The other conser- 
vation equations can be treated in a similar way. This set of equations 
could be written in the form: 


n 

a 


0 


n+1 

j 


, , n ,n+l , n .n+1 , ,n ,n+l 
+ b 0. ..,+c 0,. +d 0. 

ij i»J-l 


+ e 0 


n ^n+1 


i-lj j 


= f 


n 


(2-6) 


Here, a^ , b^ , etc., are matrices whose order is equal to the number of 
dependent variables. 0^^^^ is a column vector whose components are the 
dependent variables, f^ is a column vector containing the same number of 
elements as there are dependent variables . 

A set of multidimensional equations such as (2-6) is usually very 
time consuming to solve. Considerable time savings can be realized by splitting 
the multidimensional equation into a series of one-dimensional equations. 

9 

This was done in the present case by applying the Douglas-Gunn ADI method 
to each equation to obtain the intermediate steps. Each step involves implicit 
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solution in one of the coordinate directions. The solution of the resulting 
system of one-dimensional linear equations than only requires inversion of a 
tridiagonal matrix. 

To illustrate the ADI method, consider the conservation of mass eq, 

(2-5). For the one -dimensional equations in the x-direction (x-sweep), 

(n+1)- level quantities are evaluated at an intermediate level, denoted by *. 

The exception is that the argument of 6^ is evaluated at the n - level. 

The result is: 

(p - p )/At = - 6 (p u +pu -pu)- 8 (p V ) (2-7) 

X y 

Equation (2-7) is the x-sweep part of the conservation of mass equation. 

The y-sweep equation is formed from eq. (2-5) by evaluating (n+1) -level 
quantities at ** , except for the argument of 6^, in which (n+1) -level 
quantities are evaluated at *. Thus: 

** n. j..*n. n* nn. 

(p -p)/At=-8(pu +pu -pu) 

n Ti ti 

- 5 (o*V ^ ^ ^ 

Finally, it is convenient to subtract eq. (2-7) from eq. (2-8) and use the 
result as the y-sweep equation. This is: 

** . ** n n ** n n. 

(p - p )/At = - 6^(p V +pv -2pv) (2-9) 

Quantities denoted by the superscript * are intermediate results 
obtained from the solution of the set of one-dimensional (x-direction) equations 
of which eq. (2-7) is a member. The other equations in this set are formed 
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from*the other conservation equations in a manner similar to that of eq. (2-7). 
Likewise, each of the conservation equations is used to form y- sweep equations 
corresponding to eqs. (2-8) and (2-9). The use of eq. (2-9) instead of eq. (2-8) 

halves the storage requirement and reduces the computational effort. It has 

9 it* 2 n+1 ^n+1 

been shown that 0 is within order (At) of 0 , and so is taken as 0 

The complete set of x-sweep and y-sweep equations is given in Appendix B. 

The finite-difference forms of the conservation equations, such as 

eqs. (2-7) and (2-9), can be written in the following general form: 


c. 0. T . + b. 0, . + a. 0. ,, . = d, 

1 1-1, J 1 ij 1 1+1, J 1 


Y 0. . , + 0. . + Q! 0. = g. + n 

J 1,J-1 J IJ J 1,J+1 J J 


( 2 - 10 ) 

( 2 - 11 ) 


Here c , b , a , Y , P , and a are two-dimensional matrices. 0 is the 
dependent variable column matrix, and d , 5 , and T] are column matrices. 

The solution procedure for a sing le ti me step is as follows. First, 
eq. (2-10) is applied to successive rows (x-direction) to generate a set 
of coupled, one-dimensional equations. These implicit equations are arranged 

.X. 

. ^ 
into a block tridiagonal matrix for each row, and solved for 0 by a 

standard elimination technique. Second, eq. (2-11) is applied to successive 

** 

columns (y-direction) , Following the same procedure, 0 is computed at each 
grid point. 

The technique used to solve these equations is the LTJ decomposition and 
back-substitution (LUBS) method of Isaacson and Keller^^. A discussion of the 
application of this method to the present problem is given in reference 3. 
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3.0 BOUNDARY CONDITIONS 

The finite-difference equations have the same form for each of the 
interior grid points. These equations must be modified for cells which are 
adjacent to the boundaries of the computational field, however, to properly 
represent the physical boundary conditions. The particular boundary condi- 
tions are described here relative to the injector flow field sketch in 
Figure 2, 


3.1. Upper Wall (BC) 

Consider the cell with grid point ij adjacent to wall BC. (See Figure 3). 
Because this is an impermeable wall, 

(v)y_ = 0 (3-1) 

Also, to satisfy the no- slip condition, 

(u)y_ = 0 (3-2) 

In the x-momentum equation, the nondimens iona l shear stress at the wall is 



and Uj, is calculated from the Law of the Wall for a smooth wall: 


u_/u^ = (1/H) 'tn [(py/M')^j rJ + 5 (3-4) 

In applying (3-4), care was taken that (py/lJ>). . u^ R- > 30. It was also 
assumed that eq. (3-4) was valid between (y-) and (y+) . Then, 9u/9y = u^My. 
This provided a method of treating the shear stress terms in the y -momentum 


equation. Thus, 
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[<^T + + (puj)^ (3-5) 


and similarly for 
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The nondimensional pressure at the wall is equal to (P®)y • This 
quantity was evaluated by a linear extrapolation through (ij) and (ijj+1), 
so that. 


(pe) = 1.5 (pe). . - 0.5 (pe), ... 
y ij ijj+1 


(3-6) 


Because of eq. (3-2), it was also required that. 


(Su/3x) = 0 

y- 


(3-7) 


Normal derivatives at a surface were represented by a second-order 

5 

accurate, one-sided difference approximation : 


(30/3y)^_= (-8 0y_ + 9 0^^- Ay 


(3-8) 


This gives : 


(dv/3y)y_ = (9 v^^ 


(3-9) 


Because turbulent velocity fluctuations go to zero at a solid boundary, 
it was required that. 




(3-10) 


and. 


(k) = 0 

y- 


(3-11) 
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Following the approach used in eq. (3-5), the condition, 

(^u/5y)J = (pu^)^ (3-12) 

x+ 

was Imposed, and similarly for the value at x- . 

For the present calculations, it was assumed that the wall was 
adiabatic. Then, 


(qy)y_ = - (Se/3y)y_ = 0 (3-13) 

It was assumed that there is no variation in composition normal to the 
wall at the wall. Thus, 

(ac^/3y)y_ = 0 (3-14) 

The Law of the Wall was also used to determine certain quantities at 
(ij). Thus, 


(Su/5y) _ = uJk y^^ 

(3-15) 

1 2,_0.5 

"ij “ 

(3-16) 

®lj ' 

(3-17) 

3.2 Back Wall (CD) 


Wow consider the cell with grid point (ij) whose 

(x-) edge coincides 

with the back wall. (See Figure 4). 


For an impermeable wall. 


(u) = 0 

X- 

(3-18) 


and for no slip at the wall. 
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(v)^_ = 0 (3-19) 

Also from eqs. (3-18) and (3-19): 

(3u/3y)^_ = 0 (3-20) 

(5v/By)^_ = 0 (3-21) 

Using eq. (3-8) for normal derivatives at the wall gives: 

(5u/3x) = (9 u. . - u, .)/3 Ax (3-22) 

'x- ' rj 1+1, J 

and 

(Bv/Sx)^_ = (9 v^ - v^^^ )/3 Ax (3-23) 


The pressure at the back wall was obtained in the same way as on the upper 
wall by eq. (3-6). 


(pe)^. = 1.5(pe),j ^0.5(pe).^^^, 


(3-24) 


For an adiabatic wall. 




(3-25) 


Because of eqs. (3-18) and (3-19), the shear-work terms are 


zero: 


(s. u ) = 0 

• IX IX- 


(3-26) 


At the wall, the turbulent velocity fluctuations ate zero. Thus 


and. 



(3-27) 


a*i>x- 


= 0 


(3-28) 
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Also, similar to eq. (3-14), the assumption was made that. 


(3C /3x) = 0 

a X- 


(3-20) 


3.3 Centerline (DE) 

Because DE is a plane of symmetry, there is no mass flux across- it. 
Figure 5 shows a typical cell adjacent to DE. 



(3-30) 


This also gives. 


and. 


(3v/Sx)y_ = 0 (3-31) 

(3v/Sy) = (9v - V )/3Ay 


Also because of symmetry, 

(90/Sy)y_ = 0 (3-32) 

Combining eq. (3-32) with eq. (3-8) gives a formula for any dependent 
variable except v at (y-) : 

0^ = (9/8)0 - (1/8)0 (3-33) 

y- ij 

Equations (3-30) through (3-32) also give; 


s, u. ) = 0 

xy 1 y 


(3-34) 


3.4 Outflow Boundary (EF) 

For this case the x-f- cell boundary is located along EF. The approach 
followed here was the same as that recommended by Roache (reference 7, pp. 
279-281). Properties at (ij) are determined by a linear extrapolation 
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through and (i-2,j), giving: 


0 . . = 20 , ^ - 0 ^ _ . 
ij 1-1, J 1-2, j 


(3-35) 


3.5 Upper Boundary (AF) 

The cell with grid point (ij) has its y+ edge along AF, Flow 
properties at (ij) were determined by considering the straight left-running 
characteristic that passes through (ij). Properties are approximately con- 
stant along this characteristic, and can be determined by suitable inter- 
polation from the properties at (i-l,j-l) and (i,j-l), or (i-ljj-1) and 
(i-l,j). The procedure is described by Roache (reference 7, pp. 282-283). 
Referring to Figure 6a, when tan (p.^ + 8)^ 1 j 1 ^ Ay/Ax, then: 


1 ^-1 + ('f'/Ax)(0. , -0. , . J ( 3 - 36 ) 

p 1—1,1 1 Ij 1 1 1 1?1 1 


By this approximation, 0^ = ^ij' 


w = tan 
P 


90° 



Now, -w' = (Ax--t)/Ay by geometry. An expression for w^ can also be 
obtained from eq. (3-36). Equating these two and solving for L gives. 


i - (ta/Ay - 





when tan (p^ + 


> Lylhyi. 


(3-37) 



Figure 6a. 


Upper boundary AF, tan (iJi^ 1 j ] 



Figure 6b. Upper boundary AF. tan ^ ^ ^ 
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Next, consider Figure 6b, which depicts the case when 
tan + ®^i-l j-1 ^y/Ax. Let, 


w = tan (u +9) 

p p 


For constant properties along a left-running characteristic. 


0 = 0 . 1 . 1 + 




(3-38), 


Also, by geometry, w^ = (Ay-'t')/Ax. As before, these relationships yield a 
formula for -t: 


t = (iy/to - 

for tan (p.,,, + 9) < Ay/Ax. 

M 1-1, j-1 

Equations (3-37) and (3-39) are then used in eqs. (3-36) and (3-38), 
respectively, to compute the remaining properties at p. 

3.6 Inflow Boundary (AB) 

Flow properties at the inflow boundaries were specified and kept 
fixed as the solution developed in time. 
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4.0 RESULTS 

In this section some results of computations using the present numeri- 
cal method are given. Two sets of calculations are described. The first is 

for a supersonic, laminar base flow at conditions previously computed by 
5 

Allen and Cheng . The second is for a supersonic, turbulent base flow at 

12 

the conditions of the experiments of Lewis and Chapkis. 

The laminar flow case is simpler, as it involves fewer differential 
equations (four) . These laminar flow calculations have been very important 
to the program development, because they have served to establish effects of 
variations in boundary conditions, grid size, time step, and other parameters, 
without the additional complications of turbulent flow. Some of this work is 
still continuing, as indicated in what follows. These calculations have also 
served to establish that the numerical method is stable and convergent. 

The laminar flow calculations were made for = 3, ~ 550, and 

nondimensional initial boundary- layer thickness 6 = 0.41. The initial 
u-profile was the pol 3 momial in y used by Allen and Cheng, The correspond- 
ing V was obtained from the continuity and x-momentum boundary- layer equations, 
e was given by the Busemann integral for an adiabatic wall, and p was uniform 
and equal to its freestream value. The computational mesh size was Ax = 1/6 
and Ay = 1/12. This is coarser than that used by Allen and Cheng, and was 
chosen to increase economy while retaining reasonable accuracy during program 
development. The computational time step used was four times the Courant- 
Friedricks-Lewy (CFL) time step, which is the maximum time step for stable 
calculation using an explicit, finite-difference method. Stable and convergent 
solutions were obtained with the larger time step using the present implicit 
method, thus demonstrating potentially large savings in computer time. 

Figures 7-13 show computed results for the laminar flow case 



21 


described in the preceding. Figure 7 shows the field of velocity vectors. 

All the main features of the flow are evident here, including the initial 
boundary layer, the expansion and turning at the corner, the weak recompression 
shock and turning downstream, the recirculation region near the base, and 
the retarded flow near the plane of symmetry. The location of' the rear 
stagnation point is about 20% closer to the base than that calculated by 

5 

Allen and Cheng. The dividing streamline separates just below the comer, 

5 

in agreement with previous experiments and calculations. 

The streamline plot of Figure 8 shows the flow turning due to expan- 
sion and recompression even more clearly. Also, the recirculation region 
and the dividing streamline are indicated. The pressure contours of Figure 9 
very effectively illustrate the corner expansion and the recompression shock. 
Also note that downstream of the recirculation region there is a growing 
region where Sp/3y — 0 near the plane of symmetry. Results of this kind 
indicate the region of validity of the boundary-layer equations, which are 
usually used to describe the wake at some distance from the body. The density 
contours of Figure 10 also show the expansion and recompression. For the 
case of an adiabatic wall, the density near the wall is less than the free- 
stream value, and this is indicated in the solution. Corresponding internal 
energy contours are shown in Figure 11. 

The computed centerline pressure distribution is given in Figure 12. 

5 

Also shown are the calculated results of Allen and Cheng. The two are in 
close agreement. 

Centerline Mach number distribution is shown in Figure 13. This 
figure illustrates the extent of the recirculation region and the accelera- 
tion of the flow in the downstream direction. Downstream of the region of 
the present computation, the centerline Mach number will exceed one. A short 
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distance away from the centerline at the outflow boundary, the flow is already 
supersonic. It is perhaps worth noting that there is no singularity en- 
countered, such as the Crocco-Lees singularity, associated with the Mach one 
condition. 

Information concerning the convergence and the computation time of 
the present method is given in Figure 14. To evaluate the degree of convergence 
of the solution, the fractional change in each variable between successive 
time steps is evaluated at each grid point. The maximum fractional change for 
each variable and its location are then printed for each time step. Experience 
with the present method indicates that the density variation is a sensitive 
measure of convergence. Of course, by examining the changes in all the 
variables, as is being done, larger variations in some other variable can be 
easily detected. Figure 14 shows the fractional change in density plotted 
against number of time steps and against total computation time. Up to 
step 110, the maximum density change shown occurred on the centerline at the 
downstream boundary. Starting with step 111, the maximum change occurred on 
the back wall at about 0.6 H from the centerline. The maximum change remained 
close to the back wall for the rest of the computation. At step 450 the 
computation was stopped when the maximum fractional change in density had 
fallen to 0.0001. 

The computational time step was four times the CFL time step. This 

value was chosen arbitrarily, as being significantly larger than the explicit 

stability limit, yet not so large as to cause excessive problems related to 

9 10 

the choice of initial conditions. Briley and McDonald ’ used much larger 
time steps, but confined their attention to a much simpler flow. It is 
believed that the present resultsprovic^e a significant demonstration of the 
basic stability and convergence of the method. Other tests, including the 
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use of widely varying initial conditions and different time step sizes, must 
be performed in order to develop a more complete picture of the behavior of 
the method. This work is still being performed. The present computations 

i 

had initial p and e equal to freestream values everywhere except in the 
inflow boundary layer. For y > H, u = 1 (freestream velocity) except in 
the inflow boundary layer, and for y < H , u = 0.3. Finally, v = 0 
everywhere initially, except in the inflow boundary layer. 

The net mass flux through the field of computation was also monitored 
in these calculations. In the steady state, the net mass flux should be zero. 
At time step 450, where the calculations were halted, the total inflow differed 
from the total outflow by 1.87o. 

It seems significant that while conditions near the expansion corner 
converged most slowly, the remainder of the flow was almost unchanged after 
about time step 100. In a complex flow such as this, with a number of 
different length scales and velocity scales, different parts of the flow will 
converge to the steady state at different rates. Thus, it may be possible 
to shorten the calculation in cases where the most slowly converging portion 
is not of primary interest. 

The potential advantages in computation time of using even larger time 
steps hSVe not yet been explored. This is now being done, as reduced compu- 
tation time is probably the most attractive feature of the present numerical 
method. In this connection, the second abscissa of Figure 14 is of interest. 
Here the computation (CPU) time is indicated. For laminar flow, the method 
takes about 7 milliseconds per grid point per time step on a CDG-C5fBER 74 
machine. The present grid has 1056 points. These computation times seem 
promising, especially if significantly larger time steps can be taken. 

A calculation of a supersonic, two-dimensional, turbulent base flow 
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was made at the conditions of the experiments of Lewis and Chapkis. These 

were = 4,0, = 162,000, and 6 = 0.31, The initial u-profile in the 

boundary layer at AB was a power-law profile with exponent 1/3.7. v was 

zero, e was given by the Busemann integral for an adiabatic wall,^ The 

pressure was uniform and equal to the freestream value. As before. Ax = 1/6 

and Ay = 1/12. The time step was again four times the CFL value. In this 

2 

case, however, the k and e equations were solved to give 

In the experiments^^, the expansion corner was proceeded by a 6° 
half-angle wedge. This precompression was omitted from the present calcula- 
tions for simplicity, so that the inflow would be parallel to x. Thus, 
the TL loss in stagnation pressure across the bow shock was not simulated. 

Figure 15 shows the field of mean velocity vectors for this turbulent 
flow case. All the qualitative features of the flow are present, including 
the initial boundary layer, the turning by the corner expansion, the turning 
by the recompression shock, the recirculation near the base, and the retarded 
viscous wake. The streamlines are shown in Figure 16. Again the flow 
turning because of expansion, re compress ion, and recirculation is evident. 

The streamline spacing increases during expansion and decreases during com- 
pression. The extent of the recirculation region is too small, however, 
compared with the experimental data. Also, the dividing streamline is not 
clearly shown. This may be partly due to the pressure wiggles that appear 
in Figures 17 and 22. 

Figures 17, 18, and 19 show the contours of pressure, density, and 
internal energy, respectively. These all show the correct qualitative 
behavior, in that the expansion, recompression shock, and recirculation 
are clearly shown. Figures 20 and 21 show the k and e contours, respec- 
tively. The numerical values on the figures are one hundred times the 
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Figure 15. Velocity vectors. Turbulent flow 
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Streamlines . 


Turbulent flow. 
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Figure 17 


Pressure contours. Turbulent flow. 
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DENSITY CONTOURS 

Figure 18. Density contours. Turbulent flow. 



0.655 0.999 1-343 1.686 2.030 2.374 2.717 3.061 3.405 3.748 



INTERNRL ENERGY CONTOURS 


Figure 19. Internal energy contours. Turbulent flow. 
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Figure 20. Turbulence kinetic energy contours. 
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Figure 21. Turbulence dissipation rate contours. 
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computed values, k and e take on their largest values where the shear 
stresses are largest. 

The centerline pressure distribution is shown in Figure 22, along 

12 

with the data of Lewis and Chapkis. The calculation is fairly consistent 
with the experimental data, except that it does not show the sharp rise in 
pressure near x/H = 2. This result, plus the smaller computed recirculation 
region, may be caused by the turbulence model. The effect of the turbulence 
model on mean flow properties such as these would seem to be a useful area 
for further research. 

Figure 23 shows the corresponding Mach number distribution. 

Figure 24 is a plot of ] Ap/p]^^^ versus time step, showing the 
approach of the solution to convergence. The flow field was virtually unchanged 
after step 150. Between steps 156 and 254, j Ap/p j occurred close to the 
back wall, about halfway up. From steps 255 to 300, the location shifted 
among four places in or near the recirculation region. The calculation was 
stopped at step 300 where j Ap/p| = 0.0002. Also, at this step the net 
mass inflow differed from the net mass outflow by 2.3%. 

Careful examination of the calculations makes it clear that the pressure 
wiggles originate at the outflow boundary near the recompression shock. 

Similar behavior was previously encountered with the laminar flow calculations, 
when the condition 30 /3x = 0 was imposed along EF. This condition, when 
combined with either a linear or a quadratic variation of 0 (x) , produced 
wiggles that had their largest amplitude near the back wall. Use of equation 
(3-35) eliminated the wiggles for the laminar flow case. Further study of the 
outflow boundary condition is needed, in particular to determine whether 
the turbulence model has any effect on the appearance of wiggles. 
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Figure 23. Centerline Mach number distribution. Turbulent flow. 
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Figure 24. Variatioiv of maximum fractional change in density. Turbulent flow 
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Earlier turbulent flow calculations were made at the experimental 

6 

conditions of Badrinarayanan. -These were = 2.07, = 129,000, and 

5 = 0.55. The calculated results showed a much larger rise in pressure in the 
down-stream direction than did the experimental data. Because of the value of 
and the grid dimensions, the recompression shock passed through the upper 
boundary AF near F, instead of through EF as in the preceding cases. 

After studying the results, it became apparent that the higher pressures were 
associated with the simple-wave boundary conditions along AF. The left- 
running compression waves near the shock intersected near AF, in violation of 
the simple-wave condition. Numerical compression waves were reflected toward 
the centerline, producing an artificial pressure rise. A possible solution 
for this is to move AF upward, so that the shock exits through EF. The 
obvious disadvantage is that more grid points would be required. Alternate 
boundary conditions are also being considered. 
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5.0 CONCLUSIONS 

(1) The present numerical method computed all of the main features 

of the supersonic base flows examined, for both the laminar and turbulent cases. 
These included the comer expansion, the recompression shock, the recircula- 
tion region near the base, and the viscous wake near the centerline. 

(2) The flow calculations provided a demonstration of the basic 
stability and convergence of the numerical method. 

(3) The numerical results suggest that significant savings in 
computer time can be achieved, compared with explicit methods, by taking 
large time steps. 

(4) Further study is required to evaluate the effect of the particular 
two-equation turbulence model being used on the accuracy of the solution. In 
this connection other turbulence models should also be evaluated. 
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e 


H 

k 

K 


LIST OF SYMBOLS 

empirical coefficient in reaction-rate eq. (A20) , 
approximately 12 for the H 2 - O 2 reaction 

mean mass fraction of specie a 

mean mass fraction of fuel 

mean mass fraction of 

mean mass fraction of 0^ 

mean mass fraction of reaction product 

empirical coefficient in ^<1. (A9) , approximately 0.09 

mean mass fraction of reaction invariant 

empirical coefficients in e eq. (A16) , approximately 
1.44 and 1.92, respectively 

molecular diffusion coefficient of specie a 

mean internal energy per unit mass, made nondimens ional by e^^ 
heat of reaction per unit mass of fuel, made nondimensional by e^ 
half -height of base 

mean turbulent kinetic energy per unit mass, made nondimensional 

by 
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M 


1 


P 


Pr 




s. . 


t 


U, V 


u. 

1 


'I* 


w 

a 


x,y,z 


X, 

X 

X+ 3 X- 

y+ > y- 


Mach number, u^/ 


0.5 


mean static pressure, made nondimen sional by p, 


Prandtl number 


5e/9x. 

X 


Reynolds number, 


see eq. (A7) 

time, made nondimens ional by H/Uj^ 

mean velocity in x- and y-directions , respectively, made 
nondimensional by u^ 

mean velocity vector, made nondimensional by u^^ 
friction velocity, made nondimensional by u^ 


mean reaction rate of specie a, made nondimensional by 


Cartesian coordinates, made nondimensional by H 


x,y,z for i = 1,2,3, respectively 


right and left edges, respectively, of cell containing grid point ij 


upper and lower edges, respectively, of cell containing grid point ij 
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Greek 

y 

6 . . 

3-J 

6 

X 

6 

y 

ht 


ratio of specific heats 
Krone cker delta 




'Vhj " = (0^ - 0^.)/Ay 


computational time step, made nondimensional by H/u^^ 


Ax, Ay, 


e 


K 

M' 

P 


cr 

c 


Q 

e 


Az dimensions of cell containing grid point ij, 
made nondimensional by H 

dissipation rate of turbulent kinetic energy per unit mass, 

3 

made nondimensional by u^/H 
von Karman constant, approximately 0.43 
molecular viscosity, made nondimensional by 
eddy viscosity, made nondimensional by 
mean density, made nondimensional by 
turbulent Schmidt number, approximately 0.7 
turbulent Prandtl • number , approximately 0.7 
empirical coefficient in k-eq., approximately 1.0 
empirical coefficient in e- eq. , approximately 1.3 


0 


any dependent variable 
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Subscript 

1 denotes free stream conditions 


Superscripts 

n known time level 


n+1 unknown time level 


intermediate time level after x- sweep 

time level after y-sweep, taken to be the same as (n+1) 
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APPENDIX A. Finite Difference Forms of the Conservation Equations 
Conservation of Mass 

The differential equation for conservation of mass is, using Cartesian 
tensor notation: 

Sp/St + 3pu^/Sx^ = 0 (Al) 

This equation is then integrated over a control volume (CV) which is a cell 
of dimensions Ax, Ay, Az. Thus, 

•J^ (Sp/St) dV = - JJJ (Spu^/Bx^) dV 

CV CV 

= - f pu.n.dA (A2) 

J- 11 

CS 


using Gauss' theorem. For two-dimensional flow, the area integral over the 
control surface (GS) becomes: 


T rt* rr 

l pu.n.dA = j pu dy dz + 1 pu(-l) dydz 


CS 


x+ 


X- 


+ Jj" pv dx dz + fj pv(-l) dxdz 

y+ 


(A3) 


y- 


Next, the following approximations are made: 


JJf (3p/dt)dV — (3p/St)^j .Ax.Ay.Az 
CV 

pr 

! pu dy dz = (pu) , .Ay.Az (A4) 

J i) x+ 

xh 


etc. 
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Equation (A2) then becomes: 

(Sp/at)^^ = - (1/Ax) [(pu)^ - (p«) J - (iMy)[](pv)y4. - 


or 


ap/3t = - 6^(pu) - 6y(pv) 


(A5b) 


Conservation of Momentum : 

The differential equation is; 

B(pu.)/3t = - (a/Sy ) {pu.u. + (l/'YM^)pe 6. , - (u + p,/R )s_ , 

1 J V X J 1 ij T 1 ij 

+ (2/3)pk 6^^} (A6) 

Here, 

s. . = Su./Sx. + Bu./Sx. - (2/3)8. . 9u /Sx (A7) 

ij ij jx'^xjmm 

and the Reynolds stress tensor is modeled by: 


'^T,ij ®ij ^ij 


(A8) 


where 

pk^/e (A9) 

Equation (A6) is next integrated over the control volume that is the cell 
surrounding grid point (ij). After applying Guass' theorem and making approxi- 
mations like those in eqs. (A4) , the resulting equations for x- and y-momentum 


are: 
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X -momentum: 


Bpu/3t 


6^ -[pu^ + pe/YM^ - (M'^ + (2/3) pkj- 

6^ {puv - + p,/Rp s^} 


(AlO) 


V -momentum: 


9pv/3t 


{p“v - + H/R^) s^} 

{pv^ + pe/VM^ - (p.^ + p/Rj^) + (2/3) pkj 

(All) 


From eq . (A7 ) , 


s = (4/3) Bu/3x - (2/3) 3v/3y 



= Bu/Sy + 3v/Sx 
= (4/3) 3v/Sy - (2/3) 3 u/3x 


Conservation of EnerRy 

The differential equation is: 


(3/at) [pe + p + (K/2) pu^uj 

= - (3/aXj) {pu^ [ye + + (K/2) u^uj + + p/Pr..R^)Y q^ 

- K(p^ + p/Rj^) s^^ u^ + (2/3)K p k Ujj 


+ Kpe Kp^ s. . au,/Sx. + (2/3) Kpk au./dx. 

Tijxj '^<1.1 


(A12) 



Equation (A12) can be put into finite-difference form by following 
the same sequence of steps used for conservation of mass and momentum. The 
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result is as follows: 

(b/Bt) [pe + pC^ h^ + (K/2)p(u^ + v^)] 

= - 6^ |pe [ye + h^ + (K/2)(u^ + v^)]+ + M./Pr.Rj^)yq^ 

- (2/3) Kpku^ 

■ ^f \ + v^)] + + M-/Pr.R^)yq^ 

- K(li^ + pi/R^) s^^ u^ + (2/3) K pkv]- 

+ Kpe - K(i^ s^^ + (2/3) Kpk (Su/Sx + Sv/^y) 

In equation (A13), the following have been used: 


q = - Se/3x 

X 


q^ = - 5e/Sy 


s. u, 

IX X 


s . 

ly 


u. 

1 


- |^(4/3) 3u/3x - (2/3) 3v/^yJ u + (3u/Sy + 9v/&x)v 
= T(4/3) 3v/5y - (2/3) Su/SxJ v + (Su/5y + 3v/3x) u 


(A13) 
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S. . ' Su, /Sx, 

i: 1 j 


(4/3) TcSu/Sx)^ - (bu/3x) (3v/Sy) + (Sv/3y)^ ] 


+ (Su/3y + 3v/3x)‘ 


Conservation of k 

The differential equation is: 

S(pk)/St = - (S/3xj) [pujk - 0i^/<7j^)(^k/3xj)]- pe 

+ p,^ 3u^/3xj - (2/3) pk (4-14) 

After following the sequence of steps described in the preceding, the 
corresponding finite-difference form is: 

3(pk)/3t = - (^k/3x)j 

- [pvk -“(iI^/aj^)(3k/3y)J - pc 

+ |J.„ s. . 3u./3x. - (2/3) pk (3u/3x + 3v/3y) 

^ 3 (A15) 


Conservation of e 

The differential equation is: 


3(pe)/3t = - (3/3xj)j^pu^e - (p-^/o-^) (3s/9x^)J 


+ C 


cl 


(e/k) 3u^/3xj - (2/3) pk 


2 ,, 
pc /k 


(A16) 
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The corresponding finite-difference equation is: 

a(p€)/3t = - [pue - (Se/3x)J - 6^ j^pve - (^e/Sy)] 


+ (e/k) |^[i^ s^j - (2/3) pk (Su/3x + Sv/Sy)J + pe^/k 


Conservation of Chemical Specie m 

The differential equation is: 


^(PC )/^t = - (V^x )[ pu C - {\i>Jc + pD )(5C /3x )1 + p«J 
m j'-jni ic mmjJ m 


The corresponding finite-difference equation is: 


3(pC )/St = - 6 r puC - (M-_/cf + pD )(Sc /Sx)l 

m xLin ic mniJ 

- K [p'''=n - + % 


When C^ is C^ , can be represented by : 


~ ~ A(e/k) C C 
f f o 


When C is C. , then is zero, 

m <t> * f 

C^ is then computed from: 


C 

o 



+ X C^ 
f 


(A17) 


(A18) 


(A19) 


(A20) 


(A21) 
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where 


X == b W /a 
o 


and a and b are stoichiometric coefficients for the reaction: 


a (Fuel) + b (Oxygen) c (Product) 


For example, for the - 0^ reaction, X = 8. 

When the oxidizing medium is air, can be computed from: 




C)^^ + 2C ^ "1/^7 

o o o p J o 


((J +1.88 

p N 


Then, C is given by: 
P 


C = 1 - - 

p f 


c 

o 


N 


(A22) 


(A23) 


(A24) 


(A25) 
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APPENDIX B. The Form of the Conservation Equations for the Alternating 
Direction Implicit (ADI) Solution 

Conservation of Mass 

A backward time-difference approximation is applied to equation (A5) , 
and then the right-hand side is linearized about the nth time level. The 
result is: 

(p - p )/At =-o(pu +p u -pu) 

c . n n+1 , n+1 n n n, v 

- 6^(p v + p V - p V ) (Bl) 

This result is also given as equation (2-5) . 

For the x-sweep equations, (n+l)th level quantities are evaluated at 
an intermediate time level, denoted by * , except for arguments, of 6^. 

(p - p)/At = - 6 (p u + pu - pu) - 5 (pv) (B2) 

X y 

(x-sweep) 

Here and in what follows, quantities without a superscript are understood to 
be at the nth level. 

The y-sweep equation is formed from equation (Bl) by evaluating (n+l)th 
level quantities at ** , except for arguments of 6 , in which (n+1) level 

quantities are evaluated at * . 

(p “ p)/At = - 6 (p u + pu - pu) - 5 (p'” V + pv”" - pv) (B3) 

X y 

The final y-sweep equation is obtained by subtracting equation (B2) .from 
equation (B3) . 
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, ** * , , ** ** 

(p - p )/At =-6^(p v+pv - 2pv) (B4) 

(y- sweep) 

Equations (B2) - (B4) also appear as equations (2-7) - (2-9), respectively. 

Quantities with superscript ** are considered to be the same as the 
corresponding quantities at the (n+l)th level. 

Conservation of x-Momentum 

After time differencing and time linearization, equation (AlO) 
becomes: 

(p u + pu - 2 pu)/At 

= - 6^ [p’^‘*'^u^ + 2 puu^'*'^ - 2 pu^ + e + pe""*"^ - pe)/YM^^ 

- + H/Rp + (2/3) k + - pk)] 

- 6^ l^p uv + pu V + puv - 2puv - J (B5) 

Note that and JJ. are evaluated at the nth level,, and thus treated 

explicitly rather than implicitly. The x-sweep equation is; 

(p u + pu - 2 pu)/At 

f^2 ^ 2 

= -6 Tpu + 2 puu - 2pu + (pe+pe - pe)/VM.. 

X ^ ^ -L 

- (M-^ + lJ./R^)|^(4/3)6^u* - (2/3)6^v)J + (2/3)(p*k + pk* - pk) | 

“ [puv - (|J.^ + m,/Rj^) s^J (B6) 
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Mixed derivatives 6^ 6^ are treated explicitly at the nth level, 
because they are awkward to treat implicitly. The y-sweep equation is; 


-kit * * 

(p u+pu -pu-pu )/At 


r **,**,** „ 

=s-6 j^p uv+Pii ■'7 +puv - 3puv - 


(B7) 


Conservation of y-Momentum 

The time-linearized difference equation is: 

, n+1 . n+1 „ . , . 

(p V + pv - 2pv)/At 


= - 8^ ^p^^^^v + pu^"*'^v + puv^"*"^ - 2 puv - (p.^ + p/Rj^) 

- 8^ j^p^”*"^v^ + 2 pw^"*"^ - 2 pv^ + (p^^^e + pe^^^ - pe)/yM^ 

- (M<t + P-ZR^s + (2/3)(p'''^^k + pk"'‘^^ - pk)] (B8) 

j» yy 


The X“sweep equation is: 
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* * 

(p V + pv - 2 pv)/At 

- r * * ■* * “1 

~ ^ ^ ~ ^ ■'■ )J 

- 6 -[pv^ + pe/VM^ - (M-y ■+ p,/Rj^)[(4/3)6^v - (2/3) 6^u] 

+ (2/3) pkj (B9) 


The y-sweep equation is: 


** ** ■* * 

(p v+pv -pv-pv)/At 


r ** 2 „ ** 2 ** ** 2 

o ^ ^ P'^ - 3 pv + (p e + pe - 2 pe)/YM- 


- (ti^ +li/Rj^) (4/3) 8^(v - v) + (2/3) (p k + pk - 2 pk)J 


(BIO) 
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Conservation of Energy 

The time-linearized difference equation is: 

(1/At) [e + + (K/2) (u^ + v^)] + 

+ Kp(uu’^^^ + - 2p(e + hj^) - (3/2)Kp (u^ + v^)j'‘ 

= - 6^ {(p"'‘^^u + [ye + + (K/2) (u^ + v^)] 

+ pu (ye^"^^ + h + Kuu^^^ + Kvv^^^) 

f K. 

- pu [ 2 (ye + h^) + (3/2) K (u^ + v^)] 

•n-4»1 

- + P-/Pr*R^) Y 5^ e 

- K(ii^ + P-/R^) u^ + (2/3) K pku| 

- 6^ + pv^”^^) l^ye ■*■ ■*■ (K./2)(u^ + v )J 

, n+1 , _ n+1 , , ^ n+1 n+lv 

+ pv (ye ^R ^ 

- pv [2 (ye + h^) + (3/2)K (u^ + v^)] 

- P-/Pr.R^) y 

- K(p-^ + P-/R^) s^y u^ + (2/3) K pkv| 

+ Kpe - K iJt.„ s. . 3u./Sx. + (2/3) Kpk (6 -u + 8 v) (BH) 

i ij X j X y 

The production and dissipation terms have been treated explicitly at 


the nth level. 
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The x-sweep energy equation is: 

(1/At) {p* [e + + (K/2) (u^ + v^)] + p (e* + C* 

+ Kp (uu + vv ) - 2p(e + h^) - (3/2) Kp (u^ + v^)^ 

= - 8^ "[^(p u + pu ) I^Ye + + (K/2) (u^ + v^)J 

+ pu(Ye + c” h + Kuu + Kw ) 
r K 

- pu [ 2 (Ye + h^) + (3/2)K (u^ + v^)] 

- CJ-^/c^h V 6^ e* 

- K(H^ + p,/E^) s.^ u^ + (2/3) Kpkuj 
- 6^ (pv [ye + h^ + (K/2) (u^ + v^)] 

- K(p,^ + ]X/R^} s.y u^ + (2/3) Kpkvj 

+ Kpe - K}i s. , 5u,/Sx, + (2/3) Kpk (6 u + 8 v) 
Tijij "^'x y 


(B12) 
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The y-sweep energy equation is : 

(1/At) -|^(p - p ) j^e + + (K/2) (u^ + + P - e ) 

Vv* -* r ** * ** * 

+ P(C^ - C^) + Kp I u(u - u ) + v(v - V ) 

r A* it:k P 9 9 ~\ 

= - \(P V + pv ) [ye + C^h^+ (K/2) ( + v'")] 

** ** ** ** 

+ pv(ye + Kuu + Kvv ) 

il ix 

- pvj^3(ye + h^) + 2K(u^ + v^)J 

Conservation of k 

The time-linearized difference equation is: 

(p^*^^ k + pk"-'*'^'- 2pk)/At 

~ ~ + puk’^"*'^ - 2puk - 

- 6^ + pv’^'^^k + pvk’^"’’^ - 2pvk - (P-^/cTj^) 6y 

- pe + |j, s. . hu./dx. - (2/3 )pk (6 u + 6 v) 

T 11 1 1 ^ X V 



(B13) 


(B14) 
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The x-sweep k-equation is; 


* * 


(p k + pk - 2pk)/At 


uk + pu k + puk - 2puk - J 


5 jpvk - (ti„/cr, ) 5 k - pe+M.„ s.. Bu./3x. 
y L T k y J ^ T ij i j 


- (2/3) pk (6 u + 6 v) 

X y 


The y- sweep k-equation is: 

t ** * *5V 1 , 

(p - p )k + p(k - k)J/At 

= - 6 ^ j^p vk + pv k + pvk - 3 pvk - (iO,^/CT^) 6 ^(k - k)J 

Conservation of e 

The time-linearized difference equation is: 

(p’^'*'^ e + - 2pe)/At 

- - 6^ ut + pu’^'*'^ € + pue’^'*'^ - 2pue - ^ 

- 6^ ve + pv’^'*"^ e + pve’^'*'^ - 2pve - (P'^/a’^) 

+ (e/k) [p-^ s^^ Su^/Sxj - (2/3) pk (6^u + S^v)] - 0^2 pe^/k 


(B15) 


(B16) 


(B17) 
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The x-sweep e-equation is: 


* * 

(p e + pe - 2pe)/At 


c r * * * , *1 

0 p ue + pu e + pue - 2pue - )6 -e 

* T e X J 

“ [pve - (^Ay/ap6ye] 


Cg (e/k) [p- s 3u /Sx - (2/3)pk (6 u + 6 v)1 - C pe^/k (B18) 


The y-sweep e-equation is: 



e + p (e 



6 

y 


t ** ** ** 

p ve + pv e + pve 


- 3pve - (|J'^/<7g)6y 


(e' 


** 



(B19) 


Conservation of Species 


The time-linearized difference equation is : 
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, n+1 „ , „ n+1 „ ^ \ » 

(P - 2pC^)/At 


- r n+1 n+1 „ , „ n+1 „ „ 

= -o p uC+pu C+puC -2puC 

Xu. 3. 3. 3. 3. 

K '' <=a “a * f"’ - 2 pv + p D^) 6 


+ P W, 


(B20) 


The x-sweep species conservation equations is : 


(P* + P C* - 2 P C^)/At 


^x [P* u + P^" + P ^ ^ P ^ ^ P K 


3 


-6 ipv C - (p>^/<y + pD)6 C +pw 
y L a T c ^ a y aj ^ a 


(B21) 


The y-sweep species conservation equation is: 


r ■** * ** * “1 , 

! (p - p ) c + p (c - c ) /At 

U 3 a 3 .. 


a r ** O J- 4 - r. ** - 3 P’^C - (P-„/c^ + p D ) 6 (C - C 

o p vC + pv C + pvC a Tc ^a y a a 

y L. a ^ a ^ a 


** 


** 




(B22) 



